/*
 * SL__MUSCL.cpp
 *
 *  Created on: 08/02/2011
 *      Author: rogerio
 */

#include "SL__functions.h"

namespace PRS{

	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
	 * 							MUSCL node slope limiter
	 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

	Node_MUSCL::Node_MUSCL(){}
	Node_MUSCL::~Node_MUSCL(){}

	void Node_MUSCL::calculateNodeSlopeLimiters(VertexData* I, VertexData* J, double &SLII, double &SLJJ){
		SLII = SLJJ =1.0;
	}


	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
	 * 							MUSCL edge slope limiter
	 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
	Edge_MUSCL::Edge_MUSCL(){}
	Edge_MUSCL::~Edge_MUSCL(){
	}

	void Edge_MUSCL::calculateEdgeSlopeLimiters(SimulatorParameters* pSimPar, EdgeData* edata, int dim, const double &Sw_I, const double &Sw_J, double& SLII,double& SLJJ, double& DSwII, double& DSwJJ){
		double delta_Sw = Sw_J - Sw_I;
		VertexData* v_0 = edata->vertex_0;
		VertexData* v_1 = edata->vertex_1;

		double edIJ[3];
		for(int i=0; i<dim; i++){
			edIJ[i] = edata->length*edata->versor[i];
		}

		// Upwind-biased Interpolations
		DSwII = 2.*dot(v_0->grad_Sw,edIJ,dim) - delta_Sw;
		DSwJJ = 2.*dot(v_1->grad_Sw,edIJ,dim) - delta_Sw;

		double ratioI = (2.*DSwII*delta_Sw + qsi) / ( DSwII*DSwII + delta_Sw*delta_Sw + qsi);
		double ratioJ = (2.*DSwJJ*delta_Sw + qsi) / ( DSwJJ*DSwJJ + delta_Sw*delta_Sw + qsi);

		switch (pSimPar->getEdgeSlopeLimitFunc()){
			case MUSCL:{
				SLII = (ratioI + fabs(ratioI) + qsi)/(1. + fabs(ratioI) + qsi);
				SLJJ = (ratioJ + fabs(ratioJ) + qsi)/(1. + fabs(ratioJ) + qsi);
//				if (edata->vertex_0->ID==4 && edata->vertex_1->ID==31){
//					printf("\n%.5E\t%.5E, %.5E\t%.5E\t%.5E\t%.5E\t%.5E\t%.5E\t%.5E\t%.5E\n",
//					v_0->grad_Sw[0],v_0->grad_Sw[1],v_1->grad_Sw[0],v_1->grad_Sw[1],Sw_I,Sw_J,SLII,SLJJ,DSwII,DSwJJ);
//					printf("\n%.5E\t%.5E\t%.5E\t%.5E\t%.5E\t\n",ratioI,ratioJ,delta_Sw,edIJ[0],edIJ[1]);
//				}
				break;
			}
			case VAN_ALBADA:{
				SLII = max(.0,ratioI);
				SLJJ = max(.0,ratioJ);
				break;
			}
			case SUPERBEE:{
				double AuxLI = max(.0, min(1.,2.*ratioI));
				double AuxLJ = max(.0, min(1.,2.*ratioJ));
				SLII = max(AuxLI, min(2.,ratioI));
				SLJJ = max(AuxLJ, min(2.,ratioJ));
				break;
			}
			case MINMOD:{
				SLII = max(.0, min(1.,ratioI));
				SLJJ = max(.0, min(1.,ratioJ));
				break;
			}
			case OSHER:{
				SLII = max(.0, min(2.,ratioI));
				SLJJ = max(.0, min(2.,ratioJ));
				break;
			}
			default:
				throw Exception(__LINE__,__FILE__,"ERROR: unknown edge limiter.\n");
		}
	}
}
